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MATHEMATICAL  FLUID  DYNAMICS  OF  PLASMA  FLOW  CONTROL 

OVER  HIGH  SPEED  WINGS 

AFOSR  CONTRACT  NO.  FA9550-07-C-0039 
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Report  ofSDBD  and  flow  control  modeling  prepared  by: 

Alexander  Fedorov  and  Victor  Soloviev 

Moscow  Institute  of  Physics  and  Technology 

Moscow,  Russia 

Objectives 

Major  objectives  of  the  project  are: 

•  Develop  theoretical  models  and  computational  tools  for  SDBD  flow  control 
applications 

•  With  these  tools,  study  feasibility  of  the  delta-wing  flow  management  by  SDBD  active 
control 

The  effort  is  focused  on  the  following  tasks: 

1 .  Theoretically  and  computationally  investigate  leading-edge  separation,  vortex 
breakdown  and  other  phenomena  associated  with  flow  symmetry  breaking  and 
decreasing  lift  to  drag  ratio  for  a  delta  wing  at  high  angles  of  attack 

2.  Theoretically  and  computationally  investigate  SDBD  physics  and  develop  a 
computational  model  predicting  plasma  effects  on  aerodynamic  flows 

3.  Verify  theoretical  and  computational  tools  of  Tasks  1-2  by  comparisons  with  available 
experimental  data  as  well  as  wind-tunnel  experiments  to  be  performed  in  ITAM  under 
EOARD  sponsorship 

4.  Study  potential  for  delta-wing  flow  control  by  active  control  of  the  pre-separating 
boundary  layer  using  SDBD  along  the  wing  leading  edges 

5.  Address  scaling  issues  to  extrapolate  wind-tunnel  tests  to  full-scale  flight  conditions 

Major  accomplishments  of  2008  effort  are  summarized  in  what  follows. 
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Status  and  Major  Accomplishments 

1.  Modeling  of  surface  dielectric  barrier  discharge 

1. 1  Background 

Interest  to  the  surface  dielectric  barrier  discharge  (SDBD)  grows  permanently  due  to 
possible  advantageous  applications  to  the  boundary-layer  flow  control.  Experimental  and 
theoretical  efforts  have  been  focused  on  parametric  studies  of  SDBD-induced  jets  inside  the 
boundary  layer  [1-8].  The  discharge  models  in  these  analyses  are  either  semi-empirical  or 
based  on  simplified  numerical  simulation  of  the  discharge  evolution.  Parameters  of  SDBD- 
induced  jets  are  postulated  rather  then  predicted  by  self-consistent  modeling  of  the 
discharge  physics.  These  papers  do  not  provide  validation  of  discharge  simulations  against 
direct  experimental  measurements  of  the  discharge  parameters.  A  major  gap  in  the  plasma 
flow  control  technology  is  the  lack  of  self-consistent  physics-based  theoretical  models  and 
robust  computational  tools  that  could  provide  adequate  prediction  of  the  heat  and 
momentum  sources  induced  by  SDBD.  The  latter  are  needed  for  flow  control  applications. 

An  attempt  to  fill  this  gap  has  been  made  in  the  framework  of  physics-based  model  [9,  10] 
elaborated  under  this  project  in  2007.  A  computational  code  predicting  the  SDBD  evolution 
in  atmospheric  air  has  been  developed  for  the  cases  of  constant  positive  and  negative 
applied  voltage.  The  numerical  results  agreed  with  available  measurements  of  the  discharge 
length  and  the  surface  charge  density.  In  accord  with  experimental  data,  it  was  shown  that 
the  discharge  evolves  as  a  streamer  for  positive  electrode  polarity.  Because  of 
computational  time  restrictions  it  was  possible  to  simulate  only  one  half  of  the  discharge 
streamer  phase.  This  corresponds  to  approximately  1%  of  the  full  SDBD  cycle  including 
both  streamer  and  relaxation  phases. 

The  relaxation  phase  starts  when  the  discharge  stops  propagating  along  the  dielectric 
surface.  During  this  phase  the  discharge  plasma  decays  due  to  recombination  and  spatial 
redistribution  in  the  drift  and  diffusion  processes.  Since  the  relaxation  phase  takes  the  most 
part  of  SDBD  cycle  time  (few  microseconds  compared  to  30-50  ns  of  the  streamer  phase), 
it  produces  the  dominant  aerodynamic  effect.  Nevertheless,  to  predict  the  SDBD 
aerodynamic  forcing  both  the  streamer  phase  and  the  relaxation  phase  should  be  simulated 
properly.  In  this  connection,  the  effort  has  been  focused  on  modifications  of  the  model 
[9,10]  to  treat  the  full  SDBD  cycle  and  predict  induced  momentum  and  heat  sources 
required  for  flow  control  applications. 

1.2  Physical  model  and  governing  equations 

The  discharge  simulation  has  been  performed  in  2-D  approximation  for  the  electrode  layout 
and  the  coordinate  system  shown  in  Fig.  1.  Hereafter  the  top  semi-infinite  electrode  has  a 
height  he  =  0. 1  mm  above  the  dielectric  surface. 

The  model  accounts  for  three  types  of  charged  particles:  electrons  ( ne ),  negative  O'  ions 
(n_  )  and  positive  ions  (ni).  The  positive  ions  are  close  to  equilibrium  resulting  only 

in  one  ty  pe  of  effective  positive  ion  [9,10].  The  transport  equations  for  the  charged  particles 
are 
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(1) 


+  div J ,  =  kme  - krn  n,  -  kdnent  +  Sph ,  J,  =  n,K, E, 

r\ 

— -  +  div J_  =  0.22  katNne  —kdtn_N-kn_nt,  J_  =  -n_K  E, 
dt 


(2) 


dn. 


dt 


f-  +  divJe  =  kiNne-kdrnenl  -0.22katNne+  kdtn  N  +  Sph ,  J e  = -DeVne  -  neKeE ,  (3) 


where  J  =  flux  of  particles,  D,  K  —  diffusion  coefficient  and  mobility  respectively, 
subscripts  e,+,—  denote  electrons,  positive  and  negative  ions,  E  =  electric  field. 


electrode  <z>  = 


discharge 


electrode  (p  =  0 

Fig.  1  SDBD  electrode  layout 


The  transport  equations  are  complemented  by  Poisson  equation  for  the  if-field  potential  (p 


lS.(p  =  —4  n:e{ni  —n_—ne)  for  gas  region  (y  >  0), 
A (p  -  0  for  dielectric  layer  (- d  <  y  <  0), 

E  =  V  (p . 

Conditions  on  the  gas-dielectric  boundary  y  =  0,  x  >  0  are 


8<p 

d<p 

d<p 

dep 

—  £ 

dx 

0+ 

dx 

o-’  fy 

0+ 

-4 7zcr{x,t) , 


(4) 

(5) 

(6) 

(7) 


o- 


where  a  -  ae  +  cr+  =  surface  charge  density,  e  =  dielectric  permittivity.  Boundary 
conditions  for  the  electric  potential  are 

(p  =  0  for  y  =  -  d,  (p  =  V  for  0  <  y  <  he,  x  <  0,  (8) 


<P  =  K 


"i  i 

fx'il 

- - arctg 

— 

2  71 

{y)_ 

for  y  —>  oo,  x— >  ±  oo  , 


(9) 


where  he  =  upper  electrode  thickness,  d  =  thickness  of  the  dielectric  layer,  V  =  applied 
voltage.  The  condition  (9)  results  from  asymptotic  solution  of  equation  ts.<p  =  0  in  the 

semi-plane  y  >  0  with  the  boundary  conditions  (8)  for  d  — >  0  . 
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For  typical  SDBD  parameters  corresponding  to  the  discharge  thickness  ~  0.1  mm,  the 
electron  diffusion  flux  De\rne  is  2-3  orders  of  magnitude  smaller  than  the  drift  flux 

neKeE  .  These  fluxes  are  of  the  same  order  of  magnitude,  when  the  spatial  size  of  electron 

density  variation  is  ~  0.001  mm.  Therefore  the  electron  diffusion  can  be  neglected 
everywhere  except  a  steep  discharge  front.  For  ions,  the  diffusion  flux  is  4-5  orders  of 
magnitude  smaller  than  the  drift  flux  and  can  be  neglected  even  at  the  discharge  front. 
Hence  the  ion  transport  equations  are  hyperbolic  and  the  electron  transport  equation  is 
parabolic.  In  accord  to  the  type  of  equations  the  boundary  conditions  for  the  charged 
particle  concentrations  are 


ni  =ne  =  0,  n_  -  0  for  y  — >  co  and  y  =  0,  x  — >  ±  co, 

Pm 

J  =  -K  n  E  -  D  — =  0,  n,  =  0  . .  E,  <  0  -  n  ^  n 
*  ‘  '  y  ‘  Qy  ‘  ,  if  for>’  =  0,  x>0 

Oa  E  >  0 

,  n,  =0  -v 

£  ’  / 

nl  -  n_  -  ne  =  0  for  y  =  he ,  x  <  0. 


(10) 

(11) 

(12) 


The  condition  (10)  stands  for  the  charged-particle  densities  at  infinity.  The  condition  (11) 
stands  for  the  electron  and  ion  densities  on  the  dielectric  surface.  For  Ev  <  0 ,  when 

electrons  move  away  from  this  surface,  it  implies  that  their  hydrodynamic  flux  equals  to  the 
flux  from  the  dielectric  surface.  The  latter  is  assumed  to  be  zero  in  the  condition  (11).  For 
E  >  0 ,  when  electrons  seed  the  surface,  the  condition  of  rapid  electron  attachment  to  the 

surface  gives  ne  =  0 .  For  ions,  the  condition  =  0  is  valid  in  both  cases  because  the  ion 

diffusion  is  negligible  and  the  ion  flux  from  both  the  dielectric  surface  and  the  electrode 
surface  is  zero. 

The  condition  (12)  indicates  that  there  is  no  particle  flux  from  the  electrode  surface.  This  is 
always  valid  for  both  negative  and  positive  ions.  For  electrons,  the  condition  (12)  is  valid 
only  in  the  case  of  positive  electrode  polarity.  For  negative  electrode  polarity,  it  should  be 
replaced  by  the  condition  of  secondary  electron  emission  from  the  cathode  surface 

J ey  Ys^iy  ’  (12) 


where  ys  =  electron  secondary  emission  coefficient. 

The  balance  equations  for  the  surface  density  of  electrons,  <re,  and  ions,  cr+,  read 


dt 


=  -J 


dt 


(14) 


The  initial  conditions  are 


n,  =  ne  =  nin,  n_  =  0,  cre(x,0)  =  cr,(x,0)  =  0  at  /  =  0,  (15) 

where  nm  =  background  concentration  of  electrons  and  positive  ions. 

The  coefficients  kt,  kjr,  kr,  kat,  kdt  in  Eqs.  (l)-(3)  designate  air  ionization,  dissociative 
electron-ion  recombination,  ion-ion  recombination,  dissociative  electron  attachment  and 
detachment  rate  constants  corresponding  to  the  reactions 
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e  +  N2(02)->2e+N;(0;),  £,=0.78^+0.22^  =10 


-7.95- 


38.22 


cm  /s,  (16) 


e  +  O*  —>  0  +  0  ,  kdrl 


2x10 


-7 


300 

UK) 


,0.7 


cm3/s, 


(17a) 


6  +  04  — ^  02  +  02 , 
e  +  02  — 1 0_  +  O , 


kdrl  =  1-4x1  0“6 


f  300  N 


0.5 


-10.21-— 

r  — 3 


TM) 

cm3/s, 


cm3/s, 


(17b) 

(18) 


0“  +  N2  — >  e+  N20  ,  kdt  =  9.2-10~l3cm3/s,  (19) 

A++B"  +  M— >A  +  B  +  M  (M  =  02,N2),  £,=  2xl0'6(300/ TO1-5  cm3/s,  (20) 

where  y  —  reduced  electric  field  E/N  in  units  of  10"  V  cm  .  The  electron  temperature  Te 
and  the  ion  temperature  7)  are  measured  in  Kelvin,  and  the  Wannier  expression  is  used  to 
calculate  T,  [10].  The  appropriate  rate  constants  and  expressions  for  Ke,  De  and  Te  were 
taken  from  Refs.  [11-13].  The  ionization  rate  constants  and  the  electron  drift  velocity  were 
further  corrected  to  fit  to  the  data  [14]. 


The  source  term  Sph  in  the  right-hand  side  of  Eqs.  (1)  and  (3)  refers  to  gas  photoionization 
by  UV  radiation  from  the  discharge  region.  In  air,  the  photo-ionization  of  02  molecules  is 
caused  by  UV  radiation  of  N2(£'riu,  6 11  E*,c'4lZ*)  excited  molecules  in  the  band  98.0- 

102.5  tim.  This  source  term  was  taken  from  the  model  [15]  and  its  expression  was 
discussed  in  Refs.  [9,10]. 


In  our  previous  model  [9,10]  the  transport  equations  were  integrated  using  Particle-in-Cell 
technique  [16]  for  the  drift  motion  and  the  grid  Gauss  elimination  (sweep)  technique  for  the 
electron  diffusion.  The  ion  diffusion  was  neglected.  To  resolve  strong  nonlinearity  due  to 
the  ionization  source  in  the  right-hand  side  of  Eqs.  (1)  and  (3),  an  iteration  procedure  was 
used.  Relative  accuracy  of  iterations  was  10  ".  The  Poisson  equation  for  the  self-consistent 
electric  field  potential  was  solved  using  the  improved  Gauss-Seidel  (upper  relaxation) 
method  with  relative  accuracy  10  s.  Because  of  steep  electric  field  and  electron-ion  density 
gradients  relevant  to  the  ionization  wrave  front,  the  spatial  size  of  computational  cell  should 
be  less  than  1/500  mm.  This  grid  step  was  used  for  calculations  discussed  hereafter. 


1.3  Grid  technique  for  transport  equations 

With  the  aforementioned  numerical  approach  it  was  difficult  to  compute  the  streamer 
propagation  to  distances  >  2  mm  because  of  time  consuming  algorithm.  To  accelerate 
computations  and  increase  accuracy  of  numerical  solutions,  the  Particle-in-Cell  technique 
for  the  electron  drift  transport  solution  has  been  replaced  by  the  grid  method.  In  this 
method  a  computational  domain  is  divided  by  the  grid:  (/  =  1, ...,/),  yj  -  hyj 

(J  =  1 . J),  where  hx,  hy  are  steps  in  x-  and  y-  direction,  respectively.  The  electrical 

potential  (pt  j  and  the  electrical  charge  are  determined  in  the  grid  nodes.  The  electric  field 
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values  are  calculated  in  the  mid  points  as  ex,m/2j  =  (<pMj  -  hx , 

Eyj.j+U2=(<PiJ+l-<Pij)/hy 

The  transport  equation 

^  -  div(nKE)  + . . . ,  (21) 


is  approximated  as 

(nKEXuzj-inKEXui, 


(nKEX^n-i^EyX, 


7-1/2 


+  •••  (22) 


where  the  superscript  refers  to  a  quantity  on  the  next  time  step. 

Using  “upwind”  approximation  for  the  electron  drift  flux  components  the  electron  densities 
in  the  mid  points  are  calculated  as 


n 


1*1/2, y 


I  nMj  >  for  Ex,Mn,j  >  0 

1  \j  >  f0r  Ex,MI2,j  <  0 


(23) 


n 


i  ,y -*-1/2 


«,.;*! ,  for  >1/2  >  0 

-  f0r  E ytitj+l/2  <  0 


(24) 


The  numerical  algorithm  is  implicit,  and  the  stability  condition  gives  the  following 
restriction  on  computational  time  steps 

rx  =  mm(hx/(kEx)),  ry=xmn{hyl(JkEy)),  1/r  =  \/rx  +  l/ry .  (25) 

With  the  grid  method  computations  are  quicker  and  more  accurate.  However,  calculations 
of  streamer  evolution  to  a  distance  greater  than  few  millimeters  were  interrupted  by 
numerical  instability  that  occurred  near  the  dielectric  surface  in  the  “old”  part  of  the 
streamer  body  not  far  from  the  electrode  edge  (x  «  0.1 -0.2  mm). 

The  electron  density  profiles  for  the  streamer  body  cross-section  x  =  0. 1  mm  are  shown  in 
Fig.  2.  The  closer  the  bottom  streamer  side  to  the  dielectric  surface  the  greater  the  electron 
concentration.  This  unphysical  trend  has  been  observed  in  numerical  solutions  obtained 
with  the  Particle-in-Cell  method  as  well  as  the  grid  method.  This  could  be  caused  by 
incorrect  modeling  of  ionization  at  a  sharp  front  of  electron-density  and  electric-field 
distributions.  In  the  discussed  herein  model,  the  ionization  source  is  a  function  of  the  local 
reduced  electric  field  E/N.  If  the  streamer  bottom  side  is  close  to  the  dielectric  surface, 
then  electrons  move  against  the  Afield  force  due  to  strong  diffusion  associated  with  a  high 
concentration  gradient  and  enter  to  the  region  of  strong  £-field.  In  this  region,  the  predicted 
ionization  source  is  very  high  and  the  electron-ion  density  grows  dramatically.  The  real 
ionization  source  can  not  be  so  large,  because  the  electrons  lose  their  energy  moving 
against  the  £-field  force  and  can  not  ionize  gas  molecules  so  effectively.  Accordingly,  a 
correct  ionization  model  should  account  for  variations  of  the  electron  energy  in  the 
diffusion  process. 
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Fig.  2  Instantaneous  electron-density  njtto  profiles  at  the  streamer  cross-section  .v  =  0.1  mm  according 
to  numerical  solution  reported  in  Refs.  |9,  10];  n0=0.82xl012  cm'3. 


1.4  Refinements  of  physical  model 

1.4.1  Corrections  of  ionization  source  and  diffusion  flux 

If  the  £-field  is  strongly  inhomogeneous  in  space,  then  the  ionization  rate  constant  should 
be  obtained  from  solutions  of  the  inhomogeneous  Boltzmann  kinetic  equation. 
Nevertheless,  a  reasonable  estimate  of  this  rate  constant  can  be  made  using  the  energy' 
conservation  equation  for  electrons.  This  equation  for  the  electron  average  energy  e  reads 
[17J 


”,^-(«,»,E  +  V(B,»,))v?-e(^,E  +  V(D,B,))E  =  -^,_.in,/V-n,(r„,(26) 


where  tj  =  effective  electron  energy  loss  due  to  all  inelastic  processes,  nWe!  =  total 
electron  energy  sink  associated  with  quasi-elastic  collisions,  kt  n,=  new  non-local 

ionization  rate  constant.  The  first  two  terms  in  the  left-hand  side  of  (26)  refer  to  an  average 
electron  energy  variation  due  to  spatial  drift  and  diffusion.  The  third  term  describes  the 
electron  energy  increase  because  of  the  work  produced  by  the  E-field  force.  Assuming  that 
this  energy  is  predominantly  balanced  by  inelastic  losses,  we  get  an  approximate  equation 

e(K.n,E+V(D,n,))*  =  vk,_„n,N .  (27) 

For  the  homogeneous  case,  this  balance  is  written  as 

eKeneE2  =  rjkineN  (28) 
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and  contains  the  local  ionization  rate  constant  kt  used  for  our  previous  modeling.  From 
(27)  and  (28)  we  can  express  the  non-local  ionization  rate  constant  as 


K  „=*, 


1  + 


EVtjVv) 

KnE 2 


(29) 


With  this  correction,  numerical  solutions  become  smooth  and  stable.  The  instantaneous 
electron-concentration  and  electric-field  profiles  are  shown  in  Figs.  3  and  4.  The  streamer 
bottom  side  moves  toward  the  dielectric  surface  much  slower  than  in  our  previous 
simulations  based  on  the  local  approximation  of  ionization  rate  constant. 

The  electric-field  and  electron-concentration  gradients  have  opposite  signs  near  the 
dielectric  surface.  The  electron  diffusion  coefficient  De  increases  as  the  £-field  grows. 

Accordingly,  electrons  following  the  diffusion  flux  penetrate  to  the  region  with  larger  Dc . 


Fig.  3  Instantaneous  electron-density  profiles  at 
the  streamer  cross-section  x  =  0.5  mm,  numerical 
solution  with  the  ionization  rate  constant  (29), 
n0=0.82xlQ12  cm'3. _ 


Fig.  4  Instantaneous  />field  profiles  at  the 
streamer  cross-section  x  =  0.5  mm,  numerical 
solution  with  the  ionization  rate  constant  (29), 
E0=40.35  kV/cm. 


Sharp  £-field  gradients  and  dependence  of  the  diffusion  coefficient  on  E  could  cause 
another  computational  instability  associated  w-ith  incorrect  modeling  of  the  diffusion  flux  in 
Eq.  (3):  J djf  =  -DVne .  The  correct  diffusion  flux  is  expressed  as  [17] 

V(ZVO.  (30) 

The  new  model  comprising  the  ionization  rate  constant  (29)  and  the  diffusion  flux  (30) 
allows  us  to  simulate  the  discharge  evolution  from  the  beginning  to  the  time  instant  when 
the  discharge  stops  propagating  along  the  dielectric  surface.  This  was  demonstrated  for 
both  negative  and  positive  electrode  polarities. 

1.4.2  Correction  of  the  boundary  condition 
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The  boundary  condition  (11)  on  the  dielectric  surface  is  valid,  if  the  drift  flux  of  electrons 
or  Ey  near  the  dielectric  surface  does  not  change  sign  along  the  x-axis.  This  restriction  is 
violated  in  the  case  of  alternating  voltage  that  leads  to  a  jump  of  the  electron  concentration 
at  a  point  where  Ev  =  0  .  This  interrupts  computations,  because  the  needed  accuracy  can  not 

be  achieved  in  the  jump  vicinity.  A  general  boundary  condition,  which  is  valid  for  any 
surface,  implies  that  the  hydrodynamic  flux  equals  to  the  kinetic  flux 


-KnE. 


d(Dtnt) 

dy 


a(Ey)^~ 


+  fou 


(31) 


where 


a{Ey)  = 


(l-r)exp 


T 


for  E<  0 


(1  -r) 


f  K  E  /O  f  f  v  n  l\\ 


exp 


T 


+  4 


e  J 


1  -  exp 


K  E  A 

e  y 


T 


e  J 


(32) 


for  Ev>  0 


Vt  -  electron  thermal  velocity,  A  =  electron  mean  free  path,  r  =  surface  reflection 
coefficient  for  electrons,^,  =  electron  flux  from  the  surface.  The  latter  is  zero ,fout  =  0,  for 
the  dielectric  surface  and  the  anode  electrode,  and  it  is  fout  =  —ysJ for  the  cathode 

electrode.  The  condition  (31)  smoothly  varies  with  Ey  and,  in  the  limit  of  high  Ey,  tends  to 
the  condition  (11).  The  expression  (32)  was  derived  using  phenomenological  estimates  of 
the  electron  flux  at  a  distance  from  the  surface  being  less  than  the  mean  free  path.  These 
estimates  are  based  on  elementary  kinetic  considerations.  If  Ey  =  0  and  r  =  0,  then 

a(Ev)  =  1  and  the  electron  flux  to  the  surface  equals  to  the  thermal  flux  neVT  / 4  . 

1.5  Numerical  results  ofSDBD  modeling 

The  SDBD  evolution  has  been  computed  using  the  newr  physical  model  with  the 
aforementioned  corrections  of  the  ionization  source,  the  electron  diffusion  flux  and  the 
boundary  conditions.  The  grid  technique  described  in  Section  1.3  is  used  for  numerical 
integration  of  the  drift  transport  equations.  These  refinements  allow  us  to  simulate  both 
streamer  and  relaxation  phases  of  the  SDBD  cycle. 

1.5.1  Discharge  evolution  for  negative  electrode  polarity 

The  old  model  [9,10]  predicted  the  following  scenario  of  discharge  evolution.  First,  a  small 
cathode  region  is  formed  near  the  electrode  edge.  This  region  is  characterized  by  high 
(~1015  cm'3)  electron-ion  density  with  severe  charge  separation  caused  by  high  electric 
field.  Here  strong  ionization  takes  place.  The  bom  ions  drift  to  the  electrode  edge  whereas 
electrons  drift  to  dielectric  surface.  The  electrons  seed  the  dielectric  surface  and  create  a 
negative  surface  charge,  which  shields  the  y-component  of  external  electric  field  and  slows 
down  the  electron  motion  to  the  dielectric  surface.  In  the  nose  part  of  discharge  (near  the 
right-side  boundary  of  the  negatively  charged  part  of  dielectric  surface)  the  surface  charge 
leads  to  increasing  of  the  electric  field.  Ultimately  the  £-field  becomes  greater  than  the  air 
ionization  threshold  and  initiates  an  ionizing  wave  (streamer).  This  wave  slides  over  the 
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dielectric  surface  and  creates  a  thin  (-0.01  mm  thickness)  electron-ion  layer  on  the 
dielectric  surface  (Fig.  5a).  Accordingly,  the  old  model  predicts  that,  for  negative  electrode 
polarity,  the  discharge  evolves  as  a  near-wall  streamer.  This  contradicts  to  the  experimental 
data  [18]  showing  that  the  discharge  is  diffusive. 


6.23  ns,  (b)  - 1  =11.6  ns,  (c)  - 1  =  14.2  ns;  V-  -  4.5 
kV,  e=  8,  d—  1  mm. 


Fig.  5b  New  model  results  for  electron  density 
contours  evolution,  yellow  region  denotes  the 
electrode;  conditions  are  the  same  as  in  Fig.  5a. 


Fig.  6  Instantaneous  distributions  of  Ey(x)  on  the 
dielectric  surface:  (1)  -  t  =  0,  (2)  -  1.78  ns,  (3)  - 
3.56  ns,  (4)  -  7.12  ns;  V  =  -  4.5  kV,  e  =  8,  d  =  1 
mm,  E0  =  40.35  kV/cm. 


Fig.  7  Instantaneous  distributions  of  the  surface 
charge  density  cr(x)  in  units  of  nC/cm":  (2)  - 
1.78  ns,  (3)  -  3.56  ns,  (4)  -  7.12  ns;  V=  -  4.5  kV, 
e=  8,^=1  mm. _ 


The  new  model  predicts  that  there  is  no  ionizing  wave  moving  along  the  dielectric  surface. 
The  initial  stage  is  the  same  as  in  the  old  case  -  electrons  bom  inside  the  cathode  region 
start  to  seed  the  dielectric  surface.  However  the  electric  field  induced  in  the  discharge  nose 
is  not  high  enough  to  trigger  an  ionizing  wave  (Fig.  5b),  and  the  near-surface  layer  is 
formed  in  a  different  way:  Initially,  the  Ey  component  is  positive  on  the  dielectric  surface 
for  all  x  >  0  (curve  1  in  Fig.  6)  that  corresponds  to  zero  surface  charge.  In  accord  with  the 
boundary  condition  (31),  the  increase  of  negative  surface-charge  density  corresponding  to 
Jey< 0  continues  until  Ey  changes  its  sign.  Then  the  drift  electron  flux  becomes  positive  and 


10 


high  enough  to  compensate  the  electron  diffusion  flux  to  the  surface.  Ultimately  the  total 
flux  tends  to  zero.  Figures  6  and  7  show  temporal  evolutions  of  the  Ey  component  on  the 
dielectric  surface  and  the  surface  charge  density  at  V  =  -4.5  kV.  The  Ey  component 
becomes  negative  over  the  surface  part  with  high  enough  charge  and  ultimately  tends  to  an 
asymptotic  value  corresponding  to  zero  flux  of  electrons  to  the  dielectric  surface.  The 
instantaneous  distributions  of  Ey{x)  (Fig.  6)  correlate  with  the  surface  charge  density 
distributions  a{x)  shown  in  Fig.  7.  The  predicted  value  of  a  agrees  with  the  experimental 
data  [18]. 


In  the  aforementioned  process  of  near-surface  layer  formation,  the  electric  field  inside  this 
layer  (see  curves  3,  4  in  Fig.  6  and  Fig.  8a)  is  higher  than  the  ionization  threshold  field  E,h 
(Eth  =  32.28  kV/cm,  Eth  /E0  =  0.8  for  atmospheric  air,  Eq  -  40.35  kV/cm).  This  leads  to 
additional  ionization,  which  occurs  allover  the  layer  length  in  contrast  to  the  old  model 
predicting  ionization  in  the  streamer  head  only.  Accordingly,  the  new  model  gives  a 
diffusive  discharge  for  negative  electrode  polarity  that  agrees  with  the  experimental  data 
[18]. 


Fig.  8  SDBD  structure  in  the  near-electrode 
region  at  t  =11.6  ns:  (a)  -  EylE0  contours;  (b)  - 
electron  density  rtjn$  contours;  (c)  -  positive  ion 
density  contours;  V  =  -  5  kV,  e  =  8,  d  =  1  mm; 
«o=0.82xl012  cm'3,  E0  =  40.35  kVVcm. 


Fig.  9  SDBD  structure  in  the  discharge  nose 
region  at  t  =11.6  ns:  (a)  -  Ey/Eo  contours;  (b)  - 
electron  density'  contours;  (c)  -  positive  ion 
density  contours;  conditions  are  the  same  as  in 
Fig.  8. 


Spatial  distributions  of  SDBD  parameters  are  shown  in  Figs.  8  and  9  at  a  time  instant  that  is 
close  to  the  moment  when  discharge  stops  propagating  along  the  dielectric  surface.  These 
results  were  obtained  at  V  —  -  5  kV,  e  =  8,  d  =  1  mm,  the  grid  step  h  -  0.002  mm.  A  thin 
plasma  layer  is  formed  near  the  dielectric  surface.  The  electron  and  ion  concentrations 
inside  this  layer  approximately  equal  to  5xl0ucm'\  The  layer  thickness  smoothly 
decreases  from  0.02  mm  at  x  «  0  to  0.01  mm  at  the  discharge  nose.  There  is  no  evidence  of 
a  streamer  head  -  smooth  decreasing  of  the  electron-ion  concentration  toward  the  discharge 
nose  is  observed  (Fig.  9). 


The  cathode  layer  is  clearly  visualized  by  the  electron  density  contours  in  Fig.  8b. 
Formation  of  a  new  local  maximum  in  the  positive  ion  density  distribution  (Fig.  8c)  is 
associated  with  the  beginning  of  cathode-layer  decay.  In  this  stage,  the  electron-ion 
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recombination  inside  the  cathode  layer  becomes  greater  than  ionization,  and  the  charged 
particle  concentration  starts  decreasing.  This  process  combined  with  the  ion  drift  toward 
the  electrode  edge  results  in  the  aforementioned  local  maximum. 

Note  that  even  at  the  beginning  of  cathode-layer  decay  the  electric  field  inside  this  layer  is 
much  greater  than  the  air  threshold  ionization  value  (£)/,  =  32.28  kV/cm,  i E,h  /Et  =  0.8  for 
atmospheric  air).  Nevertheless,  the  ionization  source  is  small  because  the  electron  density 
is  low  in  the  cathode  layer. 

1.5.2  Discharge  evolution  for  positive  electrode  polarity’ 

For  positive  electrode  polarity,  the  streamer  evolution  predicted  by  the  new  model  is  also 
different  from  that  reported  in  Refs.  [9,10],  As  contrasted  to  the  old  model,  the  streamer 
does  not  ‘touch’  the  dielectric  surface  (Fig.  10).  Its  velocity  approximately  2.8  times  larger 
compared  to  the  old  case.  The  streamer  parameters  do  not  vary  notably  after  the  streamer 
formation  -  the  electron  density  contours  for  the  streamer  nose,  which  are  shown  in  Fig.  10 
at  t  —  7. 1 2  ns  and  t  =  14.2  ns,  practically  coincide. 
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Fig.  10  Streamer  evolution  predicted  b>  new  model;  electron  density  «e/«o  contours  at  V=  +4.2  kV,  e  = 
8,  d  =  1  mm. 

Although  the  streamer  does  not  touch  the  dielectric  surface,  the  surface  charge  density 
(associated  with  the  ion  drift  to  the  dielectric  surface)  decreases  only  twice  compared  to  the 
old  solution  [10].  Figure  11  shows  that  the  new  distribution  cr(x)  agrees  better  with  the 
experiment  [18]. 

The  experimentally  observed  streamer  length  (Fig.  12)  for  positive  electrode  polarity 
increases  with  the  applied  voltage.  Its  minimal  value  is  around  8  mm  under  the  near 
breakdown  condition  V  »4  kV.  Unfortunately,  Ref.  [18]  did  not  provide  the  dielectric 
layer  thickness  d  and  its  permittivity  e  for  the  data  in  Fig.  12.  In  our  calculations  we  used 
d  -  1  mm  and  e  -  8 ,  which  were  reported  in  Ref.  [18]  for  another  data  set.  In  accord  with 
experimental  observations,  the  numerical  solution  shows  that  the  streamer  breakdown  in 
atmospheric  air  has  a  threshold  nature.  If  the  applied  voltage  is  below  a  certain  level,  then 
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the  streamer  is  not  formed  and  the  discharge  just  glows  near  the  electrode  edge  and  rapidly 
decays.  According  to  our  calculations  the  breakdown  threshold  voltage  is  ~4.2  kV  for 
atmospheric  air.  At  this  voltage  the  calculated  streamer  length  is  «8  mm,  when  the 
discharge  stops  moving  along  the  dielectric  surface  and  starts  decaying.  As  shown  in  Fig. 
12,  this  prediction  agrees  well  with  the  experiment  [18]. 


Fig.  11  Surface  charge  density  a  in  units  of 
nC/cm2  for  streamer  evolution;  (1)  -  (  =  7.12  ns, 
(2)  -  t  =  14.2  ns,  (3)  -  experiment  [18);  V  -  +4.2 
kV,  e=  8,  d  =  1  mm. 


Fig.  12  Experimentally  observed  SDBD  length  in 
atmospheric  air  at  different  applied  voltage  [18]: 
(1)  -  constant  voltage  of  positive  polarity,  (2)  - 
alternating  voltage,  (3)  -  constant  voltage  of 
negative  polarity.  Symbols  show  predictions  of  the 
new  model  at  K=+4.2  kV  and  V=-5  k V  (*'=£/)• 


For  constant  applied  voltage,  the  streamer  length  is  so  large  that  it  is  extremely  time 
consuming  to  simulate  the  streamer  relaxation  phase.  To  shorten  a  streamer  length  we 
reduced  the  applied  voltage  after  the  streamer  has  been  formed  (Fig.  13). 


Fig.  13  Applied  voltage  profile  used  for  shortening  of  the  streamer  length. 

1.6  Modeling  of  relaxation  phase 

The  relaxation  phase  begins  when  the  discharge  plasma  shields  the  external  electric  field  to 
a  value  less  than  the  air  ionization  threshold  and/or  the  electron  density  becomes  small 
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inside  regions  with  high  electric  field.  In  both  cases  the  ionization  source  becomes 
negligible  in  the  discharge  region. 

It  is  prudent  to  use  different  approaches  for  simulation  of  the  discharge  formation  phase 
and  the  discharge  relaxation  phases,  since  these  phases  are  associated  with  different 
physical  processes  of  different  time  scales.  In  the  discharge  formation  phase,  the  air 
ionization  is  a  key  process  for  both  negative  and  positive  electrode  polarity.  Ionization 
creates  electron-ion  pairs  inside  a  high  electric  field  region  with  initially  available 
electrons.  The  ionized  region  polarization  in  the  external  electric  field  generates  new  spatial 
distributions  of  the  .E-field  and  the  charged  particle  concentrations.  This  process  is 
governed  by  Eqs.  (l)-(5)  without  any  simplification.  The  characteristic  time  of  the 
discharge  formation  phase  is  time  of  the  electron  drift  to  the  ionization  length  (ionization 
time).  The  typical  time  step  for  simulation  of  the  discharge  formation  is  approximately  10' 
13  s. 

In  the  relaxation  phase,  the  ionization  source  becomes  notably  less  than  the  divergence  of 
electron  drift  flux  allover  the  discharge  region 

V{neKeE)»kieffNne.  (33) 

Since  the  ionization  source  is  negligible  and  the  electron  drift  time  is  small  compared  to  the 
time  of  other  inelastic  processes  (electron-ion  recombination,  electron  detachment),  the 
transport  equation  (3)  for  electrons  splits  into  two  equations:  the  balance  equation  (34)  for 
absolute  values  of  the  electron  density  and  Eq.  (35)  governing  spatial  redistributions  of 
electrons.  Fast  electron  redistribution  occurs  in  the  electric  field  generated  by  slow  ion 
motions  governed  by  Eq.  (36).  Accordingly,  the  system  of  equations  (l)-(3)  is  reduced  to 
the  system 


=  KNne  ~  Kneni  ~ 0.22kmNnt  +  kdln_N , 
ct 

(34) 

V(DVne+neKe  E)  =  0, 

(35) 

dn 

— L  +  V («,A,E)  =  k,Nnt  -  krnenl  -  krinnt , 

Ct 

(36) 

—=-  =  0.22ka,Nne  —  k.,n  N-krin  n  . 

_  at  e  at  ri  -  i 

Ct 

(37) 

Poisson  equation  for  the  electric  potential  and  the  boundary  conditions  remain  the  same. 

The  system  (34)-(37)  is  solved  numerically  as  follows.  First,  distributions  of  the  ion 
concentration  and  corresponding  electric  field  are  calculated  with  a  time  step  relevant  to  the 
ion  drift  motion  and  managed  by  the  ion  mobility  AT,.  For  these  distributions,  Eq.  (35)  is 
integrated  using  the  relaxation  technique  with  a  much  smaller  time  step. 

Because  the  electron  mobility  Ke  is  much  larger  than  the  ion  mobility  AT/,  the  main  time  step 
for  numerical  simulation  of  the  relaxation  phase  is  Ke  /  AT  ~  500  times  greater  than  that  for 

the  discharge  formation  phase.  It  was  found  that  the  time  step  10'1"  s  is  appropriate  for 
simulations  of  the  discharge  relaxation  phase  having  the  time  scale  10'6  s. 
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1.6.1  SDBD  relaxation  for  negative  electrode  polarity 

For  V=-5  kV  the  SDBD  relaxation  process,  which  starts  with  the  charged  particle  and  the 
electric  field  distributions  relevant  to  the  end  of  discharge  formation  (Figs.  8  and  9),  is 
illustrated  by  Figs.  14-16.  The  beginning  of  relaxation  (t  =  35.6  ns)  is  shown  in  Fig.  14. 
Compared  to  the  distributions  in  Figs.  8,  9  the  electron  density  decreases  approximately  3 
times  and  the  electron  cloud  is  repulsed  from  the  dielectric  surface. 
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Fig.  14  Electron  density  njtt0  contours  at  /  =  35.6  ns  relevant  to  the  beginning  of  relaxation  phase,  V=  - 
5  kV,  s=  8,  d~  1  mm. 


Fig.  15  Density  contours  for:  (a)  -  electrons  ttjn^  (b)  -  positive  ions  and  (c)  -  negative  ions  nJn0  at 
t  =  0.623  mcs;  V  -  -  5  kV,  e  =  8,  d  =  1  mm. 

The  electron  and  ion  density  distributions  at  t  =  0.623  mcs  are  shown  in  Fig.  15  for  the 
near-electrode  region  and  the  discharge  front.  All  charged-particle  concentrations 
significantly  decrease  due  to  recombination,  and  the  electron  density  becomes  less  than 
densities  of  positive  and  negative  ions.  As  contrasted  to  the  discharge  formation  phase, 
when  plasma  was  located  inside  a  thin  (^0.02  mm)  layer  over  the  dielectric  surface,  the 
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charged  particle  cloud  moves  away  from  the  dielectric  surface  to  a  distance  « 0.1  -0.2  mm. 
As  a  result,  the  volumetric  force  acting  on  air  due  to  motions  of  charged  particles,  occupies 
much  thicker  layer  than  in  the  discharge  formation  phase.  The  volumetric  force  is  defined 
as 


f  =  e(nl  —n_  —  ne)  E , 


(38) 


and  the  power  input  per  unit  volume  is  estimated  as 

W=  j  E. 


(39) 


The  x-  andy-components  of  time-averaged  force  (in  units  of  10h  din/cm3)  are  shown  in  Fig. 
16.  Here  the  red  color  corresponds  to  the  positive  force  direction  and  the  blue  color  -  to 
negative. 
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Fig.  16  Contours  of  time-averaged  momentum  source  components  in  units  of  106  din/cm’  at  /  =  0.623 
mcs,  V=- 5  kV,  £  =  8,  d  =  1  mm;  (a)  -x-component,  (b)  -y-component. 

1.6.2  SDBD  relaxation  for  positive  electrode  polarity 


The  applied  voltage  variation  shown  in  Fig.  13  allows  us  to  slow  down  the  streamer 
propagation  and  stop  it  at  a  distance  ~  2  mm.  Further  relaxation  of  the  streamer  plasma  is 
illustrated  in  Fig.  17  at  the  time  instant  t  =  0.9  mcs.  Concentrations  of  all  charged  particles 
decreased  by  two  orders  of  magnitude  compared  to  the  streamer  formation  phase,  they 
approximately  equal  to  5xl0ncm  \  The  electron  and  positive  ion  densities  are  close  to 
each  other,  while  the  negative  ion  density  is  three  orders  of  magnitude  less.  Accordingly, 
the  negative  charged  particle  composition  of  decaying  plasma  is  opposite  to  the  case  of 
negative  electrode  polarity  at  which  the  negative  ion  density  was  greater  than  the  electron 
density. 


1  '/h 

Figures  18  and  19  show  the  time  averaged,  /  =  —  J  fjns[ (t)dt ,  momentum  and  heat  source 

tfin  0 

contours  relevant  to  the  £-field  and  charged  particle  distributions  shown  in  Fig.  17.  The  red 
color  in  Fig.  18  corresponds  to  positive  volumetric  force,  the  blue  color  -  to  negative.  The 
sign  of  the  x-component  of  momentum  source,  Fx,  is  the  same  as  in  the  case  of  negative 
electrode  polarity  because  ^-component  has  the  same  sign.  This  is  due  to  the  fact  that  the 
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potential  inside  the  streamer  body  is  approximately  equal  to  the  electrode  potential  by  the 
end  of  streamer  forming.  After  decreasing  of  the  applied  voltage  (Fig.  13)  the  electrode 
potential  becomes  less  than  the  potential  inside  the  plasma  region  (Fig.  17d)  and  the  Ex- 
component  is  negative.  Coincidence  of  the  force  direction  with  the  £-field  direction 
indicates  that  positive  ions  primarily  contribute  to  the  momentum  source. 
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Fig.  17  Contours  of:  (a)  -  electron  density  njn^  (b)  -  negative  ions  density  njn0 ,  (c)  -  positive  ions 
density  njn 0,  (d)  -  potential  (pl<p0{(pQ=  1.091  kV),  and  (e)  -  electric  field  EIEq  at  t  =  0.9  mcs  for 
alternating  voltage  shown  in  Fig.  13,  s=  8,  d—  1  mm. 
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Fig.  18  Contours  of  time-averaged  momentum  source  components  in  units  of  106  din/cnt*  at  t  =  0.9  mcs 
for  alternating  voltage  shown  in  Fig.  13,  e-  8,  d  -  1  mm;  (a)  -  .^-component,  (b)  -  ^-component. 
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Fig.  19  Contours  of  time-averaged  heat  source  in  units  of  103  \\7cm3  at  t  =  0.9  mcs  for  alternating 
voltage  shown  in  Fig.  13,  e=  8,  d=  1  mm. 
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As  contrasted  to  the  x-component  of  momentum  source,  the  ^-component,  Fy,  changes  its 
sign  when  the  electrode  polarity  is  changed.  An  absolute  value  of  both  Fx  and  Fy  is  greater 
in  the  case  of  negative  electrode  polarity.  However,  the  momentum-source  layer  thickness 
is  greater  in  the  case  of  streamer  relaxation  after  positive  applied  voltage.  Comparing  Figs. 
18  and  19  we  conclude  that  the  heat-source  layer  approximately  3  times  thinner  than  the 
momentum-source  layer. 

Evolutions  of  £-field  and  the  momentum  source  distributions  during  the  relaxation  phase  in 
the  case  of  decreasing  applied  voltage  are  shown  in  Figs.  20-22.  By  the  end  of  streamer 
forming  (/  =  17.8  ns)  and  the  beginning  of  relaxation  phase  (/  =  89  ns)  the  electric  field 
inside  the  streamer  body  is  close  to  zero.  It  starts  growing  only  after  a  notable  plasma  decay 
and  appreciable  changes  of  the  charge  particle  balance  due  to  recombination  and  drift 
motions  (see  Fig.  1 7e  for  /  =  0.9  mcs  and  Fig.  20  for  t  =  1 .6  mcs). 

The  £-field  increase  induced  by  plasma  relaxation  is  followed  by  the  momentum  source 
formation.  At  the  beginning  of  relaxation  phase  the  momentum  source  is  not  zero  only  in  a 
small  region  associated  with  the  past-streamer  head  (see  distributions  at  /  =  0.18  mcs  in 
Figs.  21  and  22).  Afterward  this  region  spreads  and  simultaneously  the  Fx  component  arises 
inside  the  streamer  body  —  the  blue  curve  appears  in  Fig.  22  at  t  =  0.54  mcs.  This  curve 
corresponds  to  Fx  =  -500  din/cm\ 

Inside  the  past-streamer  body,  the  Fx  component  is  negative,  while  in  the  region  located 
near  the  ex-streamer  head  this  force  is  positive.  For  plasma  relaxation  at  negative  polarity 
voltage,  there  is  only  a  region  with  negative  force  both  in  x-  and  y-direction.  The 
aforementioned  spatial  distributions  of  the  volumetric  force  associated  with  plasma 
relaxation  for  positive  and  negative  electrode  polarity  are  schematically  shown  by  the  blue 
dashed  regions  in  Fig.  23.  The  red  dashed  regions  show  the  heat-source  layer.  A 
combination  of  positive  and  negative  forces  shown  in  Fig.  23  may  generate  near-surface 
vortex  structures  which  could  effectively  force  the  boundary-layer  flow.  To  check  this 
assumption  the  numerically  predicted  momentum  and  heat  sources  were  approximated  by 
analytical  expressions,  which  were  used  for  calculations  of  relevant  flow  fields  (see  Section 
3).  ‘ 


x,  mm 

Fig.  20  Dynamics  of  EIEq  contours  for  a  streamer  relaxation  in  decreased  voltage  shown  in  Fig.  13;  e  = 
8,  d  =  1  mm;  yellow  dots  denote  EIE0  =  0.01  . 
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Fig.  21  Contours  of  the  >-component  of  time-  Fig.  22  Contours  of  the  ^-component  of  time- 
averaged  momentum  source  in  units  of  106  averaged  momentum  source  in  units  of  106  din/cm'; 
din/cm’  for  streamer  plasma  relaxation  in  conditions  are  the  same  as  in  Fig.  21. 
decreased  voltage;  e-  8,  d=  1  mm. 


Fig.  23  Schematics  of  momentum  and  heat  source  distributions  for  relaxation  of:  (a)  -  diffusive  plasma 
at  negative  polarity  voltage  and  (b)  -  streamer  plasma  at  positive  decreased  voltage. 
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2.  Modeling  of  the  vortex  flow  over  a  supersonic  delta  wing  at  angles  of  attack 

We  have  continued  CFD  modeling  of  the  SDBD  effect  on  aerodynamics  of  a  supersonic 
delta  wing.  The  considered  wing  configuration  is  shown  in  Fig.  24:  the  leading-edge  sweep 
angle  A  =  60° ,  the  wing  thickness  is  zero,  the  leading  and  trailing  edges  are  sharp. 
Cartesian  coordinates  (x,y,z)  are  made  nondimensional  using  the  wing  centerline  chord 

L' .  Three-dimensional  viscous  flow  past  this  delta  wing  has  been  simulated  using  3-D 
Navier-Stokes  solver  at  the  free-stream  Mach  number  A/x=1.5,  Reynolds  number 

Re,,  =  pJUlL'l pi  =  2x  10h  based  on  the  centerline  chord  L‘  =  1  m,  stagnation  temperature 
To  =  300  K  and  angles  of  attack  0°  <  a  <  30’ . 


Fig.  24  Delta  wing  configuration  and  coordinate  system. 

Numerical  solutions  of  3-D  Navier-Stokes  equations  are  performed  using  an  implicit  finite- 
volume  method.  In  all  cases  considered  hereafter,  the  flow  is  assumed  to  be  laminar.  The 
governing  equations  are  approximated  by  a  conservative  scheme.  The  flux  vector  is 
evaluated  by  an  upwind,  flux-difference  splitting  of  Roe  [19].  MUSCL  algorithm  is  applied 
with  the  third  order  TVD  space  discretization  [20].  An  Euler  implicit  discretization  in  time 
of  the  governing  equations  is  combined  with  a  Newton-type  linearization  of  the  fluxes  to 
obtain  the  system  of  algebraic  equations  [21],  This  system  is  solved  using  a  point  Gauss- 
Seidel  scheme.  The  viscosity-temperature  dependence  is  approximated  by  the  Sutherland 
law 


p'/pl  =Ty\\  +  S)/(T  +  S),  T  =  T'/T*,  5  =  110.4/7;’,  (40) 

where  asterisks  denote  dimensional  quantities,  and  temperature  T'  is  measured  in  K.  The 
fluid  is  a  perfect  gas  with  the  specific  heat  ratio  y  - 1 .4  and  Prandtl  number  Pr  =  0.72  . 

The  problem  is  solved  for  one  half  of  the  wing  with  imposing  of  the  symmetry  conditions 
on  the  plane  z  =  0 .  On  the  outflow  boundary,  the  unknown  variables  are  extrapolated  using 
the  linear  approximation.  On  the  inflow  boundaries,  the  conditions  correspond  to 
undisturbed  free  stream.  The  no-slip  boundary  conditions  are  imposed  on  the  delta  wing 
surface.  Temperature  on  the  wing  surface  equals  to  the  adiabatic  wall  temperature.  The 
computational  grid  has  approximately  4.5  xlO6  nodes.  In  the  boundary  layer  and  the  wing 
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leading-edge  region,  the  grid  nodes  are  clustered  to  increase  resolution  of  fine  flow 
structures  associated  with  boundary-layer  and  separation  flows. 

At  high  angle  of  attack  the  flow  field  contains  two  strong  vortices  generated  by  the  rollup 
of  the  shear  layer  emanating  from  the  separation  lines  located  at  the  wing  leading  edges. 
Since  the  separation  lines  are  fixed  the  major  effect  of  plasma  forcing  is  associated  with  the 
vortex  breakdown  (vortex  burst).  Detailed  numerical  study  of  the  vortex  flow  has  been 
performed  to  estimate  feasibility  of  the  vortex  breakdown  control  using  SDBD  actuators. 
The  following  configurations  have  been  treated:  the  wing-apex  SDBD,  the  leading-edge 
SDBD  and  the  multi-element  SDBD  (Fig.  25). 


Fig.  25  Red  lines  shows  SDBD  regions  for  (a)  the  wing-apex  actuator,  (b)  the  leading-edge  actuator  and 
(c)  the  multi-element  actuator. 

It  is  assumed  that  the  SDBD  actuator  generates  momentum  and  heat  sources  in  the  region 
shown  by  the  red  strip.  From  analysis  of  experimental  and  computational  studies  of  the 
SDBD  physics  Dr.  Soloviev  suggested  the  following  approximation  of  these  sources.  For 
the  wing-apex  and  multi-element  actuators,  the  volumetric  force  components  are 


(41) 


(42) 


For  the  leading-edge  actuator,  the  volumetric  force  FL  is  perpendicular  to  the  leading  edge. 
In  this  case  the  streamwise  and  spanwise  components  are  determined  as  Fx  =-Ft  cos  A, 
Fz  =  Fl  sin  A .  The  heat  source  term  in  the  energy  equation  is  approximated  as 
Q  =  (20/i  (>’)  within  each  SDBD  strip  shown  in  Fig.  25.  For  all  cases,  F0=104Ar/m3, 
y0  =  3  x  1  O'5  m ,  Q0  =  2  x  1  09  Wj m3 ,  the  width  of  SDBD  strips  is  1.5  cm. 
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2. 1  Solutions  for  the  wing-apex  actuator 

Consider  the  wing-apex  actuator  (Fig.  25a)  with  the  leading  and  trailing  edges  of  SDBD 
forcing  x,  =  3  cm  and  x2  =4.5  cm.  By  changing  sign  of  FL  we  simulate  SDBD  acting 

downstream  ( FL  >0 )  and  upstream  ( FL  <  0 ).  In  the  cases  of  a  =  0°  and  a  -  5° ,  there  is  no 

appreciable  effect  of  SDBD  on  the  flow  field  because  there  is  no  global  separation  from  the 
wing  leading  edges.  For  a  =  10°  (Fig.  26),  the  SDBD  actuator  produces  a  noticeable 
influence  on  the  wall  temperature  field  and  the  streamline  pattern.  First  evidence  of  the 
vortex  breakdown  is  observed  near  the  wing  trailing  edge  in  the  case  of  no  SDBD  forcing 
(a). 


(a)  SBD  is  off 


(b)  SBD  acts  downstream  (c)  SBD  acts  upstream 

Fig.  26  Streamlines  and  surface  temperatures  on  the  leeward  side  at  a  =  10°  . 


For  a  =  20°  (SDBD  is  off)  a  well-developed  vortex  is  observed  in  the  mid-chord  station 
(Fig.  27a).  The  burst  locus  moves  downstream  with  the  SDBD  acting  in  the  downstream 
direction  (Fig.  27b),  while  it  moves  upstream  with  the  SDBD  acting  in  the  upstream 
direction  (Fig.  27c).  This  example  demonstrates  that  the  vortex  breakdown  locus  can  be 
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controlled  by  the  SDBD  forcing  near  the  delta-wing  apex.  This  is  consistent  with  the  low- 
speed  modeling  of  Visbal  and  Gaitonde  [22],  However  the  integral  aerodynamic  forces  (lift 
and  drag  coefficients  shown  in  Figs.  28  and  29)  are  weakly  affected.  Presumably  the 
aerodynamic  loads  at  supersonic  speeds  are  not  so  sensitive  to  the  vortex  burst  locus.  For 
a  >  25° ,  the  influence  of  DBD  on  the  vortex  breakdown  is  not  so  clear,  because  the 
breakdown  point  is  very  close  to  the  wing  apex  in  all  three  cases. 

For  sufficiently  high  angles  of  attack,  the  flow  field  may  evolve  with  time  and  the 
foregoing  interpretation  should  be  taken  carefully.  To  clarify  this  issue  we  performed  direct 
numerical  simulation  of  unsteady  vortex  fields  for  the  case  of  a  -  20° .  As  expected  the 
primary  vortex  breakdown  slowly  evolves  with  time.  The  secondary  vortex  reveals  more 
unsteady  oscillatory  behavior.  This  unsteadiness  is  observed  both  with  and  without  SDBD 
forcing.  Nevertheless,  appreciable  migrations  of  the  burst  locus  wreakly  affect  the  lift 
coefficient  CL  (it  varies  in  the  range  0.836<Ci  <0.840). 


(a)  SDBD  is  off 


(b)  SDBD  acts  downstream  (c)  SDBD  acts  upstream 

Fig.  27  Streamlines  and  surface  temperatures  on  the  leeward  side  at  a  =  20°  . 
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Fig.  28  Lift  coefficient  versus  angle  of  attack. 


Fig.  29  Drag  coefficient  versus  angle  of  attack. 
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2.3  Quasi-steady  solutions  for  the  leading-edge  actuator 

Consider  the  leading-edge  SBD  actuator  shown  in  Fig.  25b.  The  SDBD  strip  coordinates 
are  specified  as  x,  =3  cm,  x2  =6  cm  and  x3  =90  cm.  Figure  30  shows  the  drag  polar,  the 

data  of  Section  2.2  for  the  wing-apex  actuator  are  also  shown  for  comparison.  In  this 
figure:  the  red  line  corresponds  to  the  case  of  SDBD  off;  the  blue  line  —  to  the  case  when 
wing-apex  SDBD  induces  the  momentum  source  in  the  positive  x-direction  (downstream 
forcing);  the  green  line  -  to  the  case  when  the  wing-apex  SDBD  induces  the  momentum 
source  in  the  negative  x-direction  (upstream  forcing);  the  magenta  square  -  to  the  case 
when  the  leading-edge  SDBD  induces  the  momentum  perpendicular  to  the  leading  edge  at 
a  =  20°.  Additional  computations  were  carried  out  to  distinguish  the  SDBD  momentum 
effect  from  the  SDBD  heating  effect.  The  black  triangle  corresponds  to  the  case  when  the 
momentum  source  is  included  into  Navier-Stokes  equations  while  the  heat  source  is  not. 
The  cyan  triangle  corresponds  to  the  opposite  situation  -  the  heat  source  is  on  while  the 
momentum  source  is  off.  In  all  cases  SDBD  forcing  produces  small  effect  on  the  wing 
aerodynamic  performance. 

Figure  31  illustrates  the  influence  of  SDBD  on  3-D  streamlines  and  the  wall  temperature  at 
the  angle  of  attack  a  =  20°.  When  SDBD  is  off  (Fig.  27a),  the  vortex  burst  is  observed  in 
the  mid-chord  station.  The  momentum  source  (Fig.  31a)  and  the  heat  source  (Fig.  31b)  lead 
to  an  appreciable  downstream  shift  of  the  vortex  burst  locus.  Nevertheless,  this  weakly 
affects  the  integral  aerodynamic  forces  (Fig.  30). 
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(a)  (b) 

Fig.  31  The  leading-edge  SDBD  effect  on  the  vortex  flow  field  at  a  =  20°  ;  (a)  the  momentum  source  is 
on  w  hile  the  heat  source  is  off;  (b)  the  heat  source  is  on  while  the  momentum  source  is  off. 


2.4  Solutions  for  the  multi-element  actuator 

Consider  the  multi-element  configuration  comprising  five  SDBD  strips  (Fig.  25c).  This 
configuration  resembles  the  SDBD  actuator  used  for  the  low-speed  wind  tunnel 
experiments  in  ITAM  [23].  The  SDBD  strip  coordinates,  xja  and  xib ,  are  given  in  Table  1. 

Figure  32  shows  the  five-strip  SDBD  effect  on  the  flow'  streamlines  and  the  w'all 
temperature  pattern.  Similar  to  the  configurations  considered  in  previous  sections,  the 
multi-element  SDBD  actuator  causes  appreciable  downstream  shift  of  the  vortex  burst. 
Nevertheless,  the  integral  aerodynamic  coefficients  vary  in  the  range  of  3%. 


Table  1 


Number  of  SDBD  strip 

x^m 

x\b>m 

i 

0.03 

0.045 

2 

0.15 

0.165 

3 

0.4 

0.415 

4 

0.6 

0.615 

5 

0.8 

0.815 
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(a)  SDBD  is  off  (b)SDBDison 


Fig.  32  The  five-strip  SDBD  effect  on  streamlines  and  surface  temperature,  a.  =  20°. 


3.  Modeling  of  SDBD  effect  on  subsonic  flow  over  a  flat  plate 

The  new  results  of  Section  1  indicate  that  SDBD-induced  momentum  and  heat  sources  have 
nontrivial  spatial  distributions,  which  are  quite  different  at  positive  and  negative  electrode 
polarities.  It  is  reasonable  to  assume  that  such  SDBD  can  produce  essentially  different 
aerodynamic  effects  on  the  near-wall  flow  depending  on  time-distributions  of  alternating 
voltage.  To  verify  this  assumption  we  have  carried  out  numerical  simulation  of  the  SDBD 
aerodynamic  performance  on  a  flat  plate  in  subsonic  free  stream. 

Consider  a  laminar  flow  over  a  flat  plate  of  length  L  =  10  cm.  The  ffee-stream  parameters 
are:  velocity  =10  m/s  ,  pressure  pM  =10'  Pa ,  density  px  =  1.2  kg/ rn  and  temperature 
Tx  =290K.  The  fluid  is  air  with  the  specific  heat  ratio  y  =  \.4  and  Prandtl  number 
Pr  =  0.72  .  The  viscosity-temperature  dependence  is  approximated  by  the  Sutherland’s  law 
(40).  The  leading-edge  of  SDBD  region  is  xx  =  2.5  cm  and  corresponds  to  ffee-stream 

Reynolds  number  ReM  =  pJJxxJpx  =  1.73x1 04  at  which  the  boundary-layer  flow  is 
laminar. 

2-D  compressible  Navier-Stokes  equations  are  solved  with  the  following  boundary 
conditions:  free-stream  conditions  on  the  inflow  and  upper  boundaries  of  computational 
domain,  ‘soft’  conditions  on  the  outflow  boundary  -  extrapolation  of  dependent  variables, 
no-slip  conditions  on  the  plate  surface,  the  wall  temperature  is  adiabatic.  The  computational 
grid  has  approximately  2.1xl05  cells.  In  the  y  -direction  there  are:  300  grid  nodes  in  the 
boundary-layer  region,  40  nodes  in  the  inviscid  flow  region.  In  the  x  -direction  there  are: 
580  grid  nodes  along  the  plate  surface,  30  grid  nodes  upstream  of  the  plate  leading  edge 
and  20  grid  nodes  downstream  of  the  plate  trailing  edge.  The  grid  nodes  are  clustered  the 
SDBD  region  2.5  cm  <  x  <  3. 1  cm  . 
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Without  SDBD  forcing  the  numerical  solution  agrees  well  with  the  Blasius  solution  for 
laminar  boundary-layer  flow  (Fig.  33). 


Fig.  33  Boundary-layer  profile  u  =U (r))  /  Ux  at  the  station  x  =  0.025  m ,  rj  =  y  /  yf^x/U^  . 


3. 1  Modeling  of  SDBD  momentum  and  heat  sources 

As  shown  in  Section  1,  the  characteristic  time-scale  of  SDBD  cycle  is  of  the  order  of  few 
microseconds.  Since  the  hydrodynamic  time-scale  is  much  longer  (~  (SDBD  length)/^  ~ 

few  milliseconds),  it  is  feasible  to  modulate  the  applied  voltage  by  a  slow  function  q>{l )  of 
hydrodynamic  time-scale  as  schematically  shown  in  Fig.  34. 


Fig.  34  Applied  voltage  consists  of  SDBD  cycles  with  time-scale  ~  few  microseconds  (blue  boxes),  which 
are  modulated  by  a  relatively  slow  function  (p(t )  (red  line)  of  hydrodynamic  time-scale  ~  few  msec. 

Using  results  of  SDBD  modeling  discussed  in  Section  1.6.2,  we  approximate  the 
volumetric  force  (Fx,Fy)  and  heat  source  Q  by  the  relations 

(*, y,  t )  =  (p{t) [Fxd  (x, y )  +  Fxh  (x, y)]  ,Fv(x,y,t)  =  (p{t)Fyd  (x, y  ) ,  (43) 
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where 


Q(x,y,t)  =  <p(t)Qd(x,t), 


(44) 


FXJ  =  AdSAx)f(y) > Fyd  =  AyjSj(x)f{y) ,  (45) 

F*  =  AxhSh{x)f{y ) ,  Qd  =  AqJgj(x)fq(y) . 


The  functions  g^x)  and  gh(x )  are  shown  in  Fig.  35,  and  the  functions  f(y)  and  f  (y)  - 
in  Fig.  36.  Their  analytical  expressions  are 
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where  ld  =  5  mm  and  lh  =  1  mm  are  lengths  of  the  streamer  body  and  the  streamer  head, 
respectively;  hd=Q.\  mm  and  hn.  =  0.005  mm  are  vertical  sizes  shown  in  Fig.  23b; 
x,  =  2.5  cm,  x2  =  x,  +  ld  and  x3  =  x2  +  lh . 


The  constants  Axd ,  A, ,  Axh  and  Aqd  in  Eq.  (45)  depend  on  the  type  of  modulation 
function  <p(t ) .  Their  values  will  be  specified  for  different  cases  considered  hereafter. 


Fig.  35  Longitudinal  distribution  of  momentum  and  heat  sources  induced  by  SDBD  at  positive  voltage 
polarity. 
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Fig.  36  Vertical  distributions  of  momentum  and  heat  sources  induced  by  SDBD  at  positive  voltage 
polarity. 

3.2  SDBD  performance  in  stationary  regime 

Consider  the  case  when  (pit )  is  a  step-function:  (p{t)  =  0  for  t  <  0 ,  (pit )  =  1  for  t  >  0 . 
During  each  SDBD  cycle  the  applied  voltage  is  positive  and  generates  a  streamer  plasma. 
The  momentum  and  heat  sources  are  distributes  as  schematically  shown  in  Fig.  23b.  After 
transient  process  associated  with  the  actuator  switching  on,  these  sources  are  averaged  over 
the  SDBD  duty  cycle  and  treated  as  steady  on  the  hydrodynamic  time  scale. 

The  constants  in  Eq.  (45)  are  evaluated  as: 

Axd=- 104  N/m  ,Ayd  =  5- 104  ;V/m3  ,Axh  =  2 -105  N/m*  ,Aqd=  2 - 1 08  If/m3 .  (51) 

The  unsteady  problem  is  solved  with  the  time  step  5  •  1 0-6  s .  The  numerical  solution  shows 
that  during  the  transient  process  a  vortex  is  generated  downstream  from  the  actuator.  The 
vortex  penetrates  outside  the  boundary  layer  and  propagates  downstream  with  the  local 
speed  of  undisturbed  flow.  Since  this  speed  quickly  increases  from  zero  to  the  free-stream 
velocity  across  the  boundary  layer,  the  vortical  structure  is  stretched  in  the  longitudinal 
direction  (Fig.  37). 

By  the  time  f  =  15x10  4  s,  the  transient  process  is  completed,  and  a  steady  flow  sets  in 
(Fig.  38).  It  is  seen  that  SDBD  actuator  performs  as  a  weakly  heated  longitudinal  jet  that 
accelerates  the  boundary-layer  flow  in  the  downstream  direction.  The  obtained  flow' 
velocity  increase  by  5  m/s  correlates  with  experimentally  observed  values  [5].  Such  a 
tangential  jet  regime  is  commonly  used  for  flow  control  applications. 


(a)  temperature  Field  (scale  from  290K  to  350K) 
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(b)  z-component  of  voticity  field  (scale  from  -10  s  s’1  to  +10  5  s'1) 


(c)  u-velocity  field  (scale  from  -1.5  m/s  to  15  m/s) 

Fig.  37  Flow  fields  at  the  time  instant  t  =  5x10  4  s  relevant  to  the  transient  process. 


(a)  temperature  field  (scale  from  290K  to  350K) 


(b)  z-component  of  voticity  field  (scale  from  -10  s  s  1  to  +10  5  s'1) 


(c)  u-velocity  field  (scale  from  -1.5  m/'s  to  15  m/s) 


Fig.  38  Steady-state  flow  fields  induced  by  SDBD,  /  =  15x10  4  s. 


3.3  SDBD  in  periodic  regime 

Consider  the  case  when  cp(t)  is  a  periodic  step-function:  ^>(/)  =  l  for  0  <t  <T  /  2 , 
(p{t)  =  -l  for  T/2<t<T,  where  the  period  is  T  =  (ld  +lh)IUx  =  6x  1CT4  s.  During  the  first 
half  of  period  =  1 ),  the  actuator  generates  a  streamer-discharge  plasma  (Fig.  23b)  and 
the  constants  in  Eq.  (45)  are  evaluated  as 

Axd  =  -3  -104  N/m\Ayd  =  104  A/m3 ,  Aa  =  2  •  104  N/m\Aqd=  107  (f/m3  .  (52) 

During  the  second  half  of  period  (<p(t)  =  — 1 ),  the  actuator  generates  a  diffusive-discharge 
plasma  (Fig.  23a)  and  the  constants  in  Eq.  (45)  are  evaluated  as 

Axd  =  -4  - 1 0s  iV/m3  ,  Ayd  =  -3  •  10s  Ar/m3  ,Axh=  0 ,  Aqd  =  108  W/m 3 .  (53) 

The  unsteady  problem  is  solved  with  the  time  step  5x10  6  s .  An  instantaneous  flow  field 
induced  by  this  actuator  is  shown  in  Fig.  39.  In  contrast  to  the  steady  regime  discussed  in 
Section  3.2,  the  SDBD  actuator  works  as  a  vortex  generator.  During  the  diffusive  discharge 
phase,  the  volumetric  force  is  relatively  larger  and  directed  upstream.  This  leads  to 
excitation  of  strong  concentrated  vortices  near  the  leading  edge  of  actuator.  These  vortices 
penetrate  to  the  outer  flow  and  propagate  downstream. 


(a)  z-component  of  voticity  field  (scale  from  -10'5  s1  to  +10  5  s'1 
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(b)  temperature  field  (scale  from  290K  to  350K) 

Fig.  39  Instantaneous  flow  field  induced  SDBD  at  t  =  20  -10  4  S  . 


Summary 
1.  SDBD  modeling 

The  physical  and  numerical  models  developed  in  2007  for  the  surface  dielectric  barrier 
discharge  (SDBD)  in  atmospheric  air  has  been  improved.  The  refined  model  and  new 
computational  code  allow  for  simulations  of  the  total  SDBD  cycle  including  the  discharge 
formation  phase  and  the  discharge  relaxation  phase. 

Numerical  studies  showed  that  the  discharge  formation  phase  lasts  few  tens  nanoseconds 
and  creates  plasma.  The  discharge  evolves  as  a  streamer  for  positive  electrode  polarity  and 
as  a  diffusive  discharge  for  negative  electrode  polarity.  The  predicted  values  of  the 
discharge  length  and  the  surface  charge  density  agree  well  with  the  available  experimental 
data  both  for  positive  and  negative  electrode  polarity.  This  validates  the  developed  physical 
model  and  computational  code.  Although  the  discharge  formation  phase  gives  negligible 
contribution  to  the  momentum  and  heat  sources  relevant  to  flow  control  applications,  it 
provides  initial  conditions  for  the  relaxation  phase. 

The  relaxation  phase  begins  when  plasma  is  redistributed  in  such  a  way  that  it  shields  the 
external  electric  field  and  reduces  the  air  ionization  rate  in  the  discharge  region  to  very 
small  value.  The  ion-electron  and  ion-ion  recombination  along  with  the  ion  drift  motion  are 
main  processes  in  the  relaxation  phase,  which  lasts  few  microseconds.  This  phase 
effectively  contributes  to  the  momentum  and  heat  sources. 

The  spatial  structure  and  charged  particle  composition  are  quite  different  in  the  discharge 
formation  and  discharge  relaxation  phases.  During  the  discharge  formation  plasma 
occupies  a  thin  near-surface  layer  of  ~0.05  mm  thickness  for  positive  electrode  polarity 
(streamer  discharge)  and  «0.02  mm  thickness  for  negative  electrode  polarity  (diffusive 
discharge).  In  the  relaxation  phase,  plasma  spreads  away  the  dielectric  surface,  and  the 
thickness  of  momentum-source  layer  increases  to  a  0.3  mm  at  both  positive  and  negative 
electrode  polarities.  For  the  case  of  streamer  relaxation  at  decreasing  applied  voltage,  the 
thickness  of  heat-source  layer  is  approximately  3  times  less  («0.1  mm),  and  for  negative 
electrode  polarity  it  is  1.5  times  less  (~0.2  mm). 
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The  momentum-source  distributions  have  a  complex  structure  in  space.  For  streamer 
relaxation  (at  positive  applied  voltage)  the  longitudinal  component  of  volumetric  force  is 
negative  in  the  ex-streamer  body  and  positive  in  the  ex-streamer  head.  The  length  of 
negative  force  region  depends  on  the  temporal  profile  and  the  applied  voltage  level  -  it 
varies  from  1  mm  to  20  mm.  The  length  of  positive  force  region  is  approximately  0.3- 1.0 
mm.  For  negative  electrode  polarity  there  is  a  region  of  negative  longitudinal  component  of 
volumetric  force  only.  The  transversal  component  could  be  positive  or  negative  inside  the 
ex-streamer  body  depending  on  the  applied  voltage  sign  and  its  temporal  profile. 

2.  Flow  control  modeling 

2. 1  Supersonic  delta  wing 

Detailed  numerical  study  of  the  vortex  flow  past  a  delta  wing  with  sharp  leading  edges  of 
60°  sweep  angle  has  been  performed  to  estimate  feasibility  of  the  vortex  flow  control  using 
SDBD  actuators.  The  following  configurations  have  been  treated:  the  wing-apex  SDBD, 
the  leading-edge  SDBD  and  the  multi-element  SDBD.  Heat  and  momentum  sources 
produced  by  these  actuators  are  modeled  by  analytical  approximations  suggested  by  Dr. 
Soloviev.  Computations  have  been  carried  out  using  3-D  Navier-Stokes  solver  for  the  free- 
stream  Mach  number  Mx  =  1.5,  Reynolds  number  Reoo=2xl06  and  angles  of  attack 

0°  <  a  <  30° . 

Since  the  boundary-layer  separation  is  fixed  (the  separation  lines  are  located  on  the  sharp 
leading  edges),  the  major  effect  of  plasma  forcing  is  associated  with  the  vortex  breakdown 
(vortex  burst).  It  was  found  that  the  vortex-burst  locus  can  be  controlled  by  the  wing-apex 
actuator.  Namely,  the  breakdown  point  moves  downstream  with  the  SDBD  acting  in  the 
downstream  direction  while  it  moves  upstream  with  the  SDBD  acting  upstream.  However, 
the  integral  aerodynamic  forces  are  weakly  affected.  The  actuator  causes  about  3% 
variations  of  the  lift  CL  and  drag  CD  coefficients.  Presumably  the  aerodynamic  loads  are 

weakly  sensitive  to  the  vortex  burst  locus  for  the  delta  wing  and  ffee-stream  parameters 
considered  herein. 

Direct  numerical  simulation  of  unsteady  flow  fields  show'ed  that  the  vortex  breakdown 
locus  evolves  with  time.  This  unsteady  behavior  is  sensitive  to  the  SDBD  forcing. 
Nevertheless,  appreciable  migrations  of  the  burst  locus  produce  small  (less  than  2%)  effect 
on  the  lift  coefficient. 

CFD  studies  showed  that  the  leading-edge  and  multi-element  SDBD  actuators  cause  an 
appreciable  downstream  shift  of  the  vortex  burst  locus.  However,  the  aerodynamic 
coefficients  are  weakly  affected. 

2.2  Subsonic  flow  over  a  flat  plate 

First-cut  numerical  simulations  of  the  SDBD  aerodynamic  performance  on  a  flat  plate  in 
subsonic  free  stream  have  been  carried  out  using  the  new  results  of  SDBD  modeling 
summarized  in  Section  1.  Numerical  solutions  of  2-D  unsteady  Navier-Stokes  equations 
showed  that  the  SDBD  actuator  can  produce  essentially  different  aerodynamic  effects  on 
the  near-wall  flow  depending  on  time-modulations  of  the  applied  voltage: 
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■  In  a  stationary  regime  with  positive  applied  voltage,  the  actuator  produces  a  weakly 
heated  longitudinal  jet  that  accelerates  the  boundary-layer  flow  in  the  dowmstream 
direction.  The  obtained  flow  velocity  increase  agrees  with  experimental  observations. 
Such  a  tangential  jet  regime  is  commonly  used  for  flow  control  applications. 

■  In  a  periodic  regime  comprising  the  diffusive  discharge  cycles  at  negative  polarity  and 
the  streamer  discharge  cycles  at  positive  polarity,  the  SDBD  actuator  predominantly 
works  as  a  vortex  generator 


3.  Future  effort 

The  obtained  results  on  SDBD  modeling  indicate  that  depending  on  temporal  profiles  of 
the  applied  voltage  it  is  feasible  to  produce  quite  different  aerodynamic  effects  on  near-wall 
flows.  These  effects  need  to  be  understood  and  sorted  out.  For  each  case  optimal  SDBD 
parameters  should  be  determined.  Scaling  issues  relevant  to  extrapolation  of  wind-tunnel 
tests  to  full-scale  flight  conditions  need  to  be  addressed.  First  this  should  be  done  for 
relatively  simple  flows,  when  the  SDBD  effect  is  ‘clean’  and  can  be  well  determined. 

In  this  connection  we  suggest  to  focus  our  next-year  effort  on  parametric  studies  of  SDBD 
performance  for  the  near-wall  flow  on  a  flat  plate.  This  basic  knowledge  will  guide  SDBD 
flow  control  strategies  for  more  complicated  aerodynamic  configurations. 

It  is  also  important  to  conduct  parametric  calculations  of  the  momentum  and  heat  sources 
induced  by  the  SDBD  actuator  and  correlate  numerical  results  by  reliable  analytical 
approximations.  The  latter  will  be  integrated  into  fluid  dynamics  tools  and  used  for  various 
flow  control  applications. 
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Personnel  Supported 

V.R.  Soloviev  and  V.M.  Krivtsov  -  SDBD  modeling 

A.V.  Fedorov  and  V.G.  Soudakov  -  flow  control  modeling 

Transitions 

Examples  of  Use  in  Technology  Applications 

•  Our  transition  prediction  techniques  are  being  applied  to  hypersonic  vehicles  by  Don 
Picetti  and  others  of  Boeing  and  are  now  being  part  of  Boeing’s  automated  design 
program  (BVIDs) 

•  Our  plasma  modeling  and  experiments  have  played  an  important  part  of  Boeing 
programs. 

Discoveries,  Inventions,  and  Patent  Disclosures 

•  “Surface  Plasma  Discharge  for  Controlling  Leading  Edge  Contamination  and  Crossflow 
Instabilities  for  Laminar  Flow,  US  6,805,325,  Bl,  issued  October  19,  2004 

•  “Surface  Discharge  for  Controlling  Vortex  Asymmetry,”  US  6,796,532  B2,  issued 
September  28,  2004 
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